March 23, 2024
This Jupyter notebook demonstrates how to define matrices in Python and then perform some basic matrix algebra.
First, we'll import numpy, which is the only package that we'll need.
import numpy as np
Here's how you can enter, for example, a 3x4 matrix. There are 3 rows of 4 elelments each (i.e. 4 columns).
m = [[1, 3, 0, 4], [2, -2, 1, 2], [-4, 1, -1, -7]]
m
[[1, 3, 0, 4], [2, -2, 1, 2], [-4, 1, -1, -7]]
The output of the matrix doesn't look very nice. We can use numpy to improve the look of the matrix output.
M = np.matrix(m)
print(M)
[[ 1 3 0 4] [ 2 -2 1 2] [-4 1 -1 -7]]
For what it's worth, I found on Stack Overflow the following only-line bit of code for formatting the matrix output in another way.
print('\n'.join(['\t'.join([str(cell) for cell in row]) for row in m]))
1 3 0 4 2 -2 1 2 -4 1 -1 -7
We can also get the shape of out matrix using numpy.
s = np.shape(M)
print(s)
print('Matrix M has', s[0], 'rows', s[1], 'columns and', np.size(m), 'elements.')
(3, 4) Matrix M has 3 rows 4 columns and 12 elements.
Another option for inputting a matrix is to first start with a $n\times m$ matrix of zeros and then assign specific values to the various elements. Below, we load a $3\times 3$ matrix of zeros.
# Make a 3x3 matrix of zeros
M1 = np.matrix(np.zeros((3,3)))
print(M1)
[[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
We can then index M1 to assign values to the various elements. Remember, Python starts counting from zero, so what we would refer to as the (1, 1) element will be accessed using M1[0, 0]. Let's make M1 a symmetric matrix.
M1[0, 0] = 1
M1[1, 1] = 4
M1[2, 2] = 8
M1[0, 1] = 3
M1[1, 0] = M1[0, 1]
M1[0, 2] = 5
M1[2, 0] = M1[0, 2]
M1[1, 2] = 7
M1[2, 1] = M1[1, 2]
print(M1)
[[1. 3. 5.] [3. 4. 7.] [5. 7. 8.]]
Of course, we can also select any element of a matrix using the same indexing method.
M1[1, 2]
7.0
We can also select entire row are columns as follows.
# Here's entire 2nd row.
row2 = M1[1]
print('row 2:', row2, '\n')
# Here's another way to get the 2nd row.
row2 = M1[1, :]
print('row 2:', row2,'\n')
# Here's the third column
col3 = M1[:, 2]
print('col 3:\n', col3)
row 2: [[3. 4. 7.]] row 2: [[3. 4. 7.]] col 3: [[5.] [7.] [8.]]
We can assign new values to an entire row.
M1[1, :] = np.arange(10, 40, 10)
print(M1)
[[ 1. 3. 5.] [10. 20. 30.] [ 5. 7. 8.]]
We can also assign new values to an entire column.
M1[:, 2] = np.transpose([np.arange(-5, -20, -5)])
print(M1)
[[ 1. 3. -5.] [ 10. 20. -10.] [ 5. 7. -15.]]
Let's now recall our $3\times 4$ matrix M and define a new $4\times 3$ matrix called N.
print('M:\n', M)
N = np.matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])
print('\nN:\n', N)
M: [[ 1 3 0 4] [ 2 -2 1 2] [-4 1 -1 -7]] N: [[ 1 2 3] [ 4 5 6] [ 7 8 9] [10 11 12]]
We can multiply our $3\times 4$ and $4\times 3$ matrices to create a new $3\times 3$ matrix. We'll use np.dot(M, N) to do the matrix multiplication $\underline{\underline{P}} = \underline{\underline{M}}\,\underline{\underline{N}}$.
P = np.dot(M, N)
print(P)
[[ 53 61 69] [ 21 24 27] [-77 -88 -99]]
We can also multiply our $4\times 3$ and $3\times 4$ matrices to create a new $4\times 4$ matrix. Numpy has another (equivalent) way to do matrix multiplication: np.matmul(M, N). Here is the result of $\underline{\underline{Q}} = \underline{\underline{N}}\,\underline{\underline{M}}$.
Q = np.matmul(N, M)
print(Q)
[[ -7 2 -1 -13] [-10 8 -1 -16] [-13 14 -1 -19] [-16 20 -1 -22]]
Note that you can also do element-by-element multiplication of two of matrices using np.multiply(). Here's an example.
x = np.matrix([[1, 2, 3], [4, 5, 6]])
y = np.matrix([[10, 20, 30], [40, 50, 60]])
z = np.multiply(x, y)
print(z)
[[ 10 40 90] [160 250 360]]
Here's some code that I found online that generates a magic square (all rows and columns add to the same number). I don't claim to understand the code.
n = 5
M5 = np.matrix([[(i+j-1+n/2)%n*n+(i+2*j-2)%n+1for j in range(n)]for i in range(n)])
print(M5)
[[11.5 13.5 20.5 27.5 4.5] [17.5 19.5 26.5 3.5 10.5] [18.5 25.5 7.5 9.5 16.5] [24.5 6.5 8.5 15.5 22.5] [ 5.5 12.5 14.5 21.5 23.5]]
Let's check if the square is indeed magic by summing all rows and columns.
for i in range(n):
print('Row', i, 'sum is', M5[i].sum())
print('Col', i, 'sum is', M5[:, i].sum())
Row 0 sum is 77.5 Col 0 sum is 77.5 Row 1 sum is 77.5 Col 1 sum is 77.5 Row 2 sum is 77.5 Col 2 sum is 77.5 Row 3 sum is 77.5 Col 3 sum is 77.5 Row 4 sum is 77.5 Col 4 sum is 77.5
Let's check the sum along the diagonal (top-left to bottom-right, called the principle diagonal). We can access the diagonal elements using np.diagonal().
dia = np.diagonal(M5)
dia_sum = sum(np.diagonal(M5))
print('The diagonal elements of M5 are:', dia)
print('The sum of the principle diagonal elements is:', dia_sum)
The diagonal elements of M5 are: [11.5 19.5 7.5 15.5 23.5] The sum of the principle diagonal elements is: 77.5
The principle diagonal has the correct sum. What about the diagonal from top-right to bottom-left, called the secondary diagonal?
We can first use np.fliplr() to reverse the order of each row. Notice that this makes the secondary diagonal of M5 become the principle diagonal of the "flipped" matrix.
M5_flipped = np.fliplr(M5)
print(M5_flipped)
[[ 4.5 27.5 20.5 13.5 11.5] [10.5 3.5 26.5 19.5 17.5] [16.5 9.5 7.5 25.5 18.5] [22.5 15.5 8.5 6.5 24.5] [23.5 21.5 14.5 12.5 5.5]]
As the result below shows, the secondary diagonal does not have the correct sum. Therefore, the matrix M5 is not a true Magic Square. When one or both of the diagonals do not have the right sum, the square is referred to as a semimagic square.
dia = np.diagonal(M5_flipped)
dia_sum = sum(np.diagonal(M5_flipped))
print('The secondary diagonal elements of M5 are:', dia)
print('The sum of the secondary diagonal elements is:', dia_sum)
The secondary diagonal elements of M5 are: [4.5 3.5 7.5 6.5 5.5] The sum of the secondary diagonal elements is: 27.5
Let's find the inverse of our $5\times 5$ semimagic square. The relevant command is np.linalg.inv().
iM5 = np.linalg.inv(M5)
print(iM5)
[[ 0.01059347 0.00226013 0.00450372 0.04104218 -0.04549628] [ 0.00290116 -0.00543218 0.05065757 -0.03588089 0.00065757] [ 0.00097808 0.04264475 -0.03780397 -0.00511166 0.01219603] [ 0.03591398 -0.03075269 0.00258065 0.00258065 0.00258065] [-0.03748346 0.00418321 -0.00703474 0.01027295 0.04296526]]
We can multiply M5 by its inverse to check if it produces the identity matrix.
# Do the matrix multiplication
I = np.dot(iM5, M5)
# Make an empty 5 x 5 matrix
Irnd = np.zeros((5, 5))
# Loop over the elements and round the values to the neareast tenth
for i in range(5):
for j in range(5):
Irnd[i, j] = round(I[i, j])
# Diplay the result in a nice format
print('\n'.join(['\t'.join([str(cell) for cell in row]) for row in Irnd]))
1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0
It's also easy to take the transpose of matrices. Of course, the transpose of a (semi)magic matrix is still (semi)magic!
M5T = np.transpose(M5)
print(M5T)
[[11.5 17.5 18.5 24.5 5.5] [13.5 19.5 25.5 6.5 12.5] [20.5 26.5 7.5 8.5 14.5] [27.5 3.5 9.5 15.5 21.5] [ 4.5 10.5 16.5 22.5 23.5]]
Someone in a quantum mechanics course may want to take the Hermitian conjugate of a matrix. First, we create a matrix with complex-valued elements.
M2 = np.matrix([[1j, 3j, 4], [2 + 3j, -2, 1 - 2j]])
print(M2)
[[ 0.+1.j 0.+3.j 4.+0.j] [ 2.+3.j -2.+0.j 1.-2.j]]
Here's the Hermitian conjugate.
M2_hc = M2.getH()
print(M2_hc)
[[ 0.-1.j 2.-3.j] [ 0.-3.j -2.-0.j] [ 4.-0.j 1.+2.j]]
Of course, we can also find the Hermitian conjugate by first evaluating the conjugate and then taking the transpose (or in the same steps in the opposite order).
M2_hc = np.conjugate(np.transpose(M2))
print(M2_hc)
[[ 0.-1.j 2.-3.j] [ 0.-3.j -2.-0.j] [ 4.-0.j 1.+2.j]]